% 此处为主函数
clc;clear;close all;
% 1. 加载活动路径
filepath = '../data/volunteer_a_all_data/csi_a1_1.dat';
% 2. 提取RSSI
% rss_sample = get_rss_from_file(filepath);
% 3. 提取CSI
csi_sample = get_csi_from_file(filepath);
% 4. 预处理
%%% 1)获取第一根天线上的数据
csi_1 = csi_sample(:,1,:);
%%% 2) 压缩维度
csi_1 = squeeze(csi_1);
%%% 3) 画图
mesh(csi_1);
%%% 4) 画图，获取30个子载波中某个子载波
csi_1_1 = csi_1(:,1);
figure;
plot(csi_1_1);title('第一个子载波对应的时域波形');
figure;
plot(csi_1);title('30个子载波对应的时域波形');
% 5.Hampel滤波、小波滤波等
x = csi_1_1; % 待处理信号
[y,i,xmedian,xsigma] = hampel(x);
n = 1:length(y);
figure;
plot(n,x);hold on
plot(n,xmedian-3*xsigma,n,xmedian+3*xsigma)
plot(find(i),x(i),'sk')
xlabel('数据包(个)');ylabel('幅度(dB)')
hold off
% legend('Original signal','Lower limit','Upper limit','Outliers')
legend('原始信号','信号下限','信号上限','异常值')
sgtitle('Hampel滤波')

len = length(y);
Fs = 320;
t0 = length(y) / Fs; % len/320(fs)
t = linspace(0,t0,len);
% 1）中值滤波比较
figure;subplot(321);
plot(y);xlabel('数据包(个)');ylabel('幅度(dB)');title('原始信号')
median_ = 5;
% way 1
% median_value = movmedian(y,median_);
% % way 2
median_value = medfilt1(y,median_); %? fail？？？
% % medfilt1，在原信号上重新绘制
% figure;plot(y);hold on; 
% plot(median_value,'r');
% sgtitle('Median Filter, k=5')
subplot(322);
plot(median_value);xlabel('数据包(个)');ylabel('幅度(dB)');title('中值滤波')
% sgtitle('Median Filter, k=5')

% 2） 均值滤波
mean_ = 5;
mean_value = movmean(y,mean_);
% % figure;
% % plot(y);hold on;plot(mean_value,'r');
% % sgtitle('Mean Filter, k=5')
subplot(323);
plot(mean_value);xlabel('数据包(个)');ylabel('幅度(dB)');title('均值滤波')
% sgtitle('Mean Filter, k=5')

% 3）平滑滤波
smooth_test = smooth(y);
subplot(324)
plot(smooth_test');xlabel('数据包(个)');ylabel('幅度(dB)');title('平滑滤波')
% sgtitle('Smooth')

% 4） butter滤波
% 设计butter低通滤波器
% [g,Wn]=buttord(Wp,ws,Rp,Rs);
fk = 1000;
f_p = 100; % 低通滤波器通带截止频率
f_s = 150; % 低通滤波器阻带截止频率 
Wp=f_p/(0.5*fk); % =20/50=0.4 ; %(通带截至频率fp = 20Hz)
Ws=f_s/(0.5*fk); % =30/50=0.6;  %(阻带截止频率fs = 30Hz)
Rp=1;  %通带最大衰减
Rs=30; %阻带最小衰减
fs = 1000; % 采样频率
[g,Wn]=buttord(Wp,Ws,Rp,Rs);
[b,a]=butter(g,Wn);
csi_size = length(y);
f=linspace(0,csi_size,csi_size); %????
k=filter(b,a,y);% 滤波
subplot(325);
plot(f,k);xlabel('数据包(个)');ylabel('幅度(dB)');title('巴特沃斯低通滤波')
set(gca,'YLim',[16 23]);


% 5) dwt 小波去噪，重构信号
%***************************************
%         wdenoise 测试
%***************************************
xden = wdenoise(y);
subplot(326)
plot(xden);xlabel('数据包(个)');ylabel('幅度(dB)');title('小波滤波')
